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Abstract 

In this work we design a new domain decomposition method for the Euler equations in 2 dimen- 
sions. The basis is the equivalence via the Smith factorization with a third order scalar equation to 
whom we can apply an algorithm inspired from the Robin-Robin preconditioner for the convection- 
diffusion equation. Afterwards we translate it into an algorithm for the initial system and prove 
that at the continuous level and for a decomposition into 2 sub-domains, it converges in 2 iterations. 
This property cannot be preserved strictly at discrete level and for arbitrary domain decompositions 
but we still have numerical results which confirm a very good stability with respect to the various 
parameters of the problem (mesh size, Mach number, . . .). 
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1 Introduction 



The need of using domain decomposition methods when solving partial differential equations is nowadays 
more and more obvious. The challenge is now the acceleration of these methods. Different possibilities 
were studied such as the use of optimized interface conditions on the artificial boundaries between subdo- 
mains or the preconditioning of a substructured system defined at the interface. The former were widely 
studied and analyzed for scalar problems such as elliptic equations in Lio90, EZ98 , for the Hclmholtz 
equation in |BD97I IUN981 IGMN01I ILVL05| convection-diffusion problems in jJNROlj . For time depen- 
dent problems and local times steps, see for instance [GHN0II |HH03 . The preconditioning methods 
have also known a wide developpement in the last decade. The Neumann-Neumann algorithms for sym- 
metric second order problems JBGLTV89 Man92l |RT91| has been the subject of numerous works, see 
TW04] and references therein. An extension of these algorithms to non-symmetric scalar problems (the 
so called Robin-Robin algorithms) has been done in [ATNVOO GGTN04] for advection-diffusion problems. 

The generalization of these kinds of methods to the system of the Euler equations is difficult in the 
subsonic case in dimensions equal or higher to two. As far as Schwarz algorithms are concerned, when 
dealing with supersonic flows, whatever the space dimension is, imposing the appropriate characteristic 
variables as interface conditions leads to a convergence of the algorithm which is optimal with regards 
to the number of subdomains. This property is generally lost for subsonic flows except for the case of 
one-dimensional problems, where the optimality is expressed by the fact that the number of iterations 
is equal to the number of subdomains (see Bj0rhus |Bj095| and Quarteroni [Qua90 for more details). 
In the subsonic case and in two or three dimensions, we can find a formulation with classical (natural) 
transmission conditions in |Qua90[ TGFS98I | QS96| or with more general interface conditions in |Gle98| and 
optimized transmission conditions in |DN04) . The analysis of such algorithms applied to systems proved 
to be very different from the scalar case, see |DLN02I IDLN 04 . As far as preconditioning methods are 
concerned, to our knowledge, no extension to the Euler equations was done. 

In this paper, we consider a preconditioning technique for the system of the compressible Euler 
equations in the subsonic case. The paper is organized as follows: in Section |2 we will first show the 
equivalence between the 2D Euler equations and a third order scalar problem, which is quite natural by 
considering a Smith factorization of this system, see WRL95 or Gan66c . In Section we define an 
optimal algorithm for the third order scalar equation. It is inspired from the idea of the Robin-Robin 
algorithm ATNVOO applied to a convection-diffusion problem. We also prove by using a Fourier analysis 
that this algorithm converges in two iterations. Afterwards in Section^ we back-transform it and define 
the corresponding algorithm applied to the Euler system. All the previous results have been obtained at 
the continuous level and for a decomposition into 2 unbounded subdomains. After a discretization in a 
bounded domain we cannot expect that these properties to be preserved exactly. Still we can show in 
Section[3]by a discrete convergence analysis that the expected results should be very good. In Sectional 
numerical results confirm the very good stability of the algorithm with respect to the various parameters 
of the problem (mesh size, Mach number, . . .). 

2 A third order scalar problem 

In this section we will show the equivalence between the linearized Euler system and a third order scalar 
equation. The motivation for this transformation is that a new algorithm is easier to design for a scalar 
equation than for a system of partial differential equations. 
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2.1 The compressible 2D Euler equations 

In the following we will concentrate ourselves on the conservative Euler equations in two-dimensions: 

dW f 
(2.1) — + V.F(W)=0 , W=(p,pV,E) t . 

at 

In the above expressions, p is the density, V = (it, v) is the velocity vector, E is the total energy per unit 
of volume and p is the pressure. In equation (|2.1I) . W — W (x, t) is the vector of conservative variables, x 
and t respectively denote the space and time variables and F(W) — (Fi(W), F2(W)) is the conservative 
flux vector whose components are given by 

Fi(W) — (pu, pu 2 + p, puv, u(E + p)) , F 2 (W) = (pv, puv, pv 2 + p, v(E +p)) ■ 

The pressure is deduced from the other variables using the state equation for a perfect gas p — (j s — 
l)(E — 5/0 || V || 2 ) where 7 S is the ratio of the specific heats (7^ = 1.4 for the air). 

2.2 Equivalence of the Euler system to a scalar equation 

The starting point of our analysis is given by the linearized form of the Euler equations (|2.1() written in 
primitive variables (p, u, v, S). In the following we suppose that the flow is isentropic, which allows us to 
drop the equation of the entropy (which is totally decoupled with respect to the others). We denote by 
W = (P, U, V) T the vector of unknowns and by A and B the jacobian matrices of the fluxes Fi(w) to 
whom we already applied the variable change from conservative to primitive variables. In the following, 
we shall denote by c the speed of the sound and we consider the linearized form (we will mark by the bar 
symbol, the constant state around which we linearize) of the Euler equations: 



(2.2) VW = — + Ad x W + Bd y W = f 

characterized by the following jacobian matrices: 
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We can re- write the system 112.211 by denoting (3 — i > under the form 
(2.4) VW = (fil + Ad x + Bd y ) W = f 

We will study this system with the help of the Smith factorization. 



2.2.1 Smith factorization 

We first recall the definition of the Smith factorization of a matrix with polynomial entries and apply it 
to systems of PDEs, see Gan66b, Gan66a, Gan98] or |WRL95j and references therein. 

Theorem 1 Let n be an integer and C an invertible n x n matrix with polynomial entries in the variable 
X : C = {cij(X))i<ij< n . 

Then, there exist three matrices with polynomial entries E, D and F with the following properties: 
• det(E)=det(F)=l 
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• D is a diagonal matrix. 

• C = EDF. 

Moreover, D is uniquely defined up to a reordering and multiplication of each entry by a constant by a 
formula defined as follows. Let 1 < k < n, 

• S k is the set of all the submatrices of order k x k extracted from C . 

• Det k = {Det{B k )\B k e S k } 

• LD k is the largest common divisor of the set of polynomials Det k . 
Then, 

LD k {\) 



(2.5) 

(by convention, LDq = 1). 



D kk {\) 



LD k _ x {\y 



\<k<n 



Application to the Euler system We first take formally the Fourier transform of the system 12. 411 
with respect to y (the dual variable is £) . We keep the partial derivatives in x since in the sequel we shall 
consider a domain decomposition with an interface whose normal is in the x direction. We note 



(2.6) 



(3 + ud x + i^v pc 2 d x 
P = | jd x P + ud x + iiv 

f + ud x + iv£ 



ipc?i 




We can perform a Smith factorization of V by considering it as a matrix with polynomials in d x 
entries. We have 



(2.7) 
where 

(2.8) 
and 



and 



(2.9) 



where 



V = EDF 



D 




E = 



(u{c 2 - u 2 )) 1 / 3 



ipc 2 £ 






u 



(3 + ud x + ivt; E 2 



F = — 



( f3 + ud x + ijv 
i^pc 2 

Ox 

pu 
u 



(3 + ud x + i^v 



pu 



i£pc 2 ) 



E 9 



P + i^v ) 

. (-uc 2 + u 3 )d xx + (2u 2 - c 2 )(f3 + itv)d x + u((P + i£,v) 2 + £, 2 c 2 ) 
' c 2 {i(3 + i^v) ' 
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(2.10) G = P + ud x + i£v 
and 

(2.11) t = f3 2 + 2i£uvd x + 2(3{ud x + i£v) + (c 2 - v 2 )£ 2 - (c 2 - u 2 )d xx 
The operators showing up in the diagonal matrix have a physical meaning: 

Q = (3 + ud x + vd y 

is a first order transport operator where the time derivative is replaced by [3 and 

C = I3 2 + 2uvd xy + 20(ud x + vd y ) - (c 2 - v 2 )d yy - (c 2 - u 2 )d xx 
is the advective wave operator where d\ is replaced by (3 for I — 1, 2. 

Equation 12.8|l suggests that the derivation of a domain decomposition method (DDM) for the third 
order operator CQ is a key ingredient for a DDM for the compressible Eulcr equations. 

3 A new algorithm applied to a scalar third order problem 

In this section we will describe a new algorithm applied to the third order operator found in section [21 
We want to solve 

(3.1) CG(Q) = g 

where Q is scalar unknown function and g is a given right hand side. The algorithm will be based on 
the Robin-Robin algorithm |AN97I [ATNV00 for the convection-diffusion problem. Then we will prove 
its convergence in 2 iterations. We first note that the elliptic operator C can also be written as: 

(3.2) C = -div(AV) + aV + f3 2 , A = ( ^ ~-f ^Ttf ) where a = ^ 

Without loss of generality we assume in the sequel that the flow is subsonic and that u > and thus we 
have < u < c. 

3.1 The algorithm for a two-domain decomposition 

We consider now a decomposition of the plane R 2 into two non-overlapping sub-domains fli — (— oo, 0) x M 
and = (0, oo, 0) x K. The interface is T = {x = 0}. The outward normal to domain fi^ is denoted rii, 
i = 1,2. Let Q l ' k , i = 1, 2 represent the approximation to the solution in subdomain i at the iteration k 
of the algorithm. We define the following algorithm: 

ALGORITHM 1 We choose the initial values Q 1 ' and Q 2 ' such that QQ 1 ' = GQ 2,a and we compute 
(Q l,k+1 )i = i.2 from (Q** )j=i 2 by the following iterative procedure: 

Correction step We compute the corrections Q 1,k and Q 2,k as solution of the homogeneous local prob- 
lems: 

CGQ 2 k = o m n 2l 

CGQ hk = in fii, 

(3.3) { { {AV - ±a)GQ 2 ' k -n 2 = 7 fc , on T, 



(AV - U)GQ l k ■ ni = 7 fc , on T. 



Q 2 k = 0, on T. 
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where i k = -\ [AVQQ 1 ^ ■ n x + A\JQQ 2 > k ■ n 2 ] . 

Update step.JVe update Q 1,k+1 and Q 2 ' k+1 by solving the local problems: 



( CgQ l k+1 = g, m O l5 
(3.4) { { GQ 2 ' k+1 =GQ 2 ' k + S k , onT, 



CGQ 2 k+1 = g, in Q 2 , 
GQ 2 ' k+1 = GQ 2 ' k + 5 k , 

Q2M+1 Q .i, . Q .i, 



where 5 = \ 



QQ 1 ^ + GQ 2 ' k ■ 
Proposition 1 Algorithm^ converges in 2 iterations. 



Proof We use the Fourier transform technique. For the sake of the analysis we consider the previous 
algorithm written in term of error vector e t,k (x,y) = (Q l,k — Q\o, i ){x,y). The error {e l ' k )i = \ : 2 satisfies 
Algorithm ^ with g = 0. 

We will first describe what happens locally inside each subdomain after proceeding to a Fourier transform 
in the y direction and then we prove the convergence of the algorithm^ by computing in Fourier space 
the effect of the correction and the update steps. We denote by e(x, £) the Fourier transform of a function 
e(x,y): 

e(x,£) = / e{x,y)e-^ y dy 



We first study solutions to the homogeneous equation CQ{ei) = in domain f2j, i = 1,2. We take its 
Fourier in the y direction and get: 

(3.5) tQe 1 = (/3 + ud x + i£v) {{u 2 - c 2 )d xx + 2u(/3 + i£v)d x + £ 2 + (/? + i£v) 2 )e l = 

We seek the solution in the form e l (x,£) ~ e A ^^.and we find three possible values for A: 



(3-6) Ai j2 (C) = , A 3 = = — • 

cr — u z u 

therefore the solution writes e l (x,£) = au{S,)e Xl ^ x + a2i(Qe X2 ^ x + aai(^)e X3 ^' x . We also impose that 
ei (resp. e 2 ) is bounded as x tends to — oo (resp. oo). Taking into account the sign of the real parts of 
(Aj)j- = i )2 ,3 it means that we have: 

(3.7) e 1 (x,0=a 1 e Xl ^ )x and e 2 (x, £) = a 2 e x ^ )x + a 3 e x ^ )x 

In order to ease the notations we call: 



(3.8) a(0 =(3 + i£v and R(£) = yj{(3 + i^v) 2 + £ 2 (c 2 - u 2 ) 

The initial guesses of the algorithm does not satisfy a specific partial differential equation. Therefore, it 
is possible to use formula (|3.7|) only for k > 1. With obvious notations, we write: 

(3.9) S 1,fc (x,0 = a k e Xl ^ x and e 2 ' fe (x,£) = a k e x ^ )x + a^ 3 ^ 
and we have also 

(3.10) e llfc (z,0 = a k e Xl ^ x and e 2 ' k (x^) = a§e Aa(fi)a + age* 8 ^* 
Using, Qe 1 ' 1 — Qe 2 ' 1 we get: 

(3.11) a r := a\(a(0c + uR(£)) = c^(a(£)c - uR(0) 
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Now we estimate the Fourier transform of the correction e. By using 1|3.1U|) in the interface conditions of 
the correction step (|3.3I) . we first get 7 1 = i£ar-R(£) and then: 



(3.12) 



_i ctr _i ar ~i ££(c — it ) ar 



a(0+uR(0' z a(Oc-UR(0' J a(£) a(Oc-ufl(0 



Now we estimate the Fourier transform of the update e. By using (|3.9|l and (|3.1U|) in the interface condi- 
tions of the update step J22J, we first get that Qe 1 ' 1 + S 1 = Qe 2 ' 1 + S 1 = and e 1 ' 1 + e 1 ' 1 = on T. This 
means that the coefficients (a 2 )i<i< 3 satisfy homogeneous systems of linearly independent equations. 
Therefore af = for i — 1,2,3 and then e 1 ' 2 = e 2,2 = 0. Therefore, the convergence is achieved in two 
steps. ■ 



4 A new algorithm applied to the Euler system 

After having found an optimal algorithm which converges in two steps for the third order model problem, 
we focus on the Euler system by translating this algorithm into an algorithm for the Euler system. It 
suffices to replace the operator CQ by the Euler system and Q by the last component F(W) 3 of F{W) 
in the boundary conditions. The algorithm reads: 

ALGORITHM 2 We choose the initial values W 1 ' and W 2 ' such that GF{W lfi ) 3 = gF(W 2 >°) 3 and 
we compute (W*' +1 )t=i,2 from (W /l ' fe )i=i.2 by the following iterative procedure: 

Correction step We compute the corrections W 1 ^ and W 2,k as solution of the homogeneous local 
problems: 



VW 2 < k = inn 



2: 



VW x ' k = m ftx, 

(4.1) { { {AV - ±a)gF(W 2 > k ) 3 -n 2 = 1 k , onT, 

(AV - \a)gF{W^ k ) 3 ■ ni = 7 fe , on Y. 

F(W 2 ' k ) 3 = 0, on T. 

where 7 fc = ~\ [AV£F(W^ fc ) 3 • ni + AVgF{W 2 > k ) 3 ■ n 2 ] . 
Update step.M^e update W 1,k+1 and W 2 ' k+1 by solving the local problems: 
(4.2) 

vw 2 - k+1 = f. m n 2 



gF{W 1 > k+1 ) 3 =gF(W 1 > k ) 3 + S k , onT. 



gp{w 2 ' k+1 ) 3 = gF{w 2 > k ) 3 + <5 fc , on r, 

F(W 2 > k+1 ) 3 = F(W^ k ) 3 + F{W l ' k ) 3 , on V. 



where S k = ■k 



gF{w lk ) 3 + gF(w 2 ' k ) 



This algorithm is quite complex since it involves second order derivatives of the unknowns in the boundary 
conditions on gF (W) 3 . It is possible to simplify it. By using the Euler equations in the subdomain, we 
have lowered the degree of the derivatives in the boundary conditions. After lengthy computations that 
we omit here, we find a simpler algorithm. We write it for a decomposition in two subdomains with an 
outflow velocity at the interface of domain but with an interface not necessarily rectilinear. In this 
way, it is possible to figure out how to use for a general domain decomposition. 

In the sequel, n = (n x ,n y ) denotes the outward normal to domain Cli, d n = V • n the normal 
derivative at the interface, d T = V • r the tangential derivative, U n = Un x + Vn y and U T = —Un y + Vn x 
are respectively the normal and tangential velocity at the interface between the subdomains. Similarly, 
we denote u n (resp. u T ) the normal (resp. tangential) component of the velocity around which we have 
linearized the equations. 



7 



ALGORITHM 3 We choose the initial values W lfi = (P 1 - , U 1 - , V 1 - ), i = 1,2 such that P 1 ' = P 2 <° 
and we compute W l,k+1 from W l,k by the iterative procedure with two steps: 

Correction step We compute the corrections W 1 ^ and W 2,k as solution of the homogeneous local 
problems: 



(4.3) 



vw lk = o, m n x , 

-09 + u T d T )U^ n + u n d T U l T - k = 7 fc , on T. 



' vw 2 ' k = 0, in Q 2 , 

{(3 + u T d T )U 2 ' k - u n d T U 2 - k = 7 fc , on T 
k P 2 ' k + pu n (j^ k = 0, on T. 



(/? + Ur0r)(E#* - U^ k ) + U n d r (U}>* - U 2 ' k ) 



where ^ k — - 

Update step. We compute the update of the solution W 1,k+1 and W 2 ' k+1 as solution of the local problems: 
(4.4) 



where S k = g 



P 



l.k 



P 



2.k 



VW 2,k+l = y 2j in 
p2,fe+l = p2.k +§ k^ m 

(P + pU n U n ) 2 > k+1 = {P + pUnUn) 1 * + {P + pUnUn) 1 *, OH T. 

2 divided into two non overlapping half planes, algorithms^ and [3] 



Proposition 2 For a domain fl = 
are equivalent and both converge in two iterations. 

Proof The Euler equations are invariant with respect to rotations so we can write at the interface the 
equations l|2.4(l in the referential given by (n, r). We simply have to replace d x (resp. d y ) by <9„(resp. d T ), 
U (resp. V") by U n (resp. U T ) and u (resp. v) by u n (resp. u T ). We denote by Q the Fourier transform 
along the interface of a function Q and £ is the dual variable. 

We first consider the boundary conditions 14.2J1 of Algorithm [3 and prove that they correspond to the 
ones (|4.4|l of Algorithm First of all, from (|2.9|l . we have: 



(4.5) F(W) 3 
We apply now the operator Q and get: 

(4.6) Q(FW)z) : 



1 



(3 + i£u T 



(u n P + pu 2 JJ n ) 



(3 + i^U T 



{Q{P) + pUnG{U n )) 



By using the Euler equations satisfied by W we can substitute G(U n ) with —l/pd n P (we can omit the 
right handside since it will appear on both sides of the boundary condition) and obtain: 



(4.7) 



G(F(W) 3 ) = - 



-(G(P)-U n d n P) = — 



(3 + i£u T j3 + i£u T 

Therefore, boundary condition of l|4.2|l on the boundary of domain ill reads: 



(iu T £ + f3)(P) = -UnP- 



(4.8) 



u n P^ k+1 = u,„ [ P 



I -It i \^p^-.k _|_ p2,k- 



By simplifying with u n , we get the boundary condition of Ij4.4[l on the pressure on the boundary of 
domain 



8 



We now show how to obtain the second boundary condition in domain both in the update step 14.4fl 
and in the correction step l|4.3|) . From (|4.5() . we infer: 

(4-9) ^—(u n P 2 > k+1 + pulTJl^ = ^_ {Un pl,k + ptfjjlk + u Jl,k + - U 2jl, k) 

p + it,u T p + iE,u T 

Multiplying by (/3 + i£u T )/i2„, we obtain the second boundary condition of l|4.3[l and of l|4.4|) . 

We now derive the first boundary condition of the correction step (|4.3|) from the corresponding boundary 

condition of (|4.1(l . We have to consider 

(4.10) B{W) := (AV - ^)GF{W) 3 
From ((S3 and lf4"?j> . we have: 

(4.11) myf) = -u n {{c 2 - u 2 n )d n - u n (P + iu T £))(P) 

In order to replace the normal derivative on P, we write the Euler system in the form: 

(4.12) d x W = -A- 1 {/3W + Bd y W-f) 

We get (once again omitting the right hand side / that will appear on both sides of the boundary 
conditions): 

(4.13) n P = _ 2 Un . [-Untf + U T d T ){P) + pd 2 (P + U T d T )(U n ) - U n pc 2 d T (U r )] 

Using this equation in (|4.1U|I , 

(4.14) B(W) = ~u n pc 2 [{l3 + u T d T )(U n ) - u n d T (U T )] 

To obtain the first boundary condition of the correction step (|4.3|l . it suffices to multiply l|4.1|l by ~u n pc 2 . 
This shows the equivalence of both algorithms. The convergence in two steps comes from the fact that 
Algorithm |21 was derived directly from Algorithm ^ which converges in two steps. M 



5 Discretization 

In this section we will first present the discretization method used, a finite volume method on a uniform 
grid. Then we propose a strategy of discretization of the boundary conditions of the algorithm [3] applied 
to the Euler system and we present some theoretical discrete estimates of the convergence of the method. 

5.1 A finite volume discretization 

We consider a domain fl and the boundary value problem associated to l|2.4|l with classical(natural) 
boundary conditions (see |DLJN04j ) on dil. This BVP is discretized using a finite volume scheme where 
the flux at the interface of the finite volume cells is computed using a Roe type solver. We recall the 
method already described in DN04 . In order to discretize the BVP we consider a regular quadrilateral 
grid where a vertex Vij is characterized by 

-J) Ax, (j-l)Ay),i,jeZ 



2 J 'V 2, 

We associate to each vertex a finite volume cell, CV,- = [(i — l)Ax , iAx] x [(j — l)Ax , jAx] which is 
a rectangle having as a center the vertex uy. A first order vertex centered finite volume formulation 
simply writes: 

1 1 1 eedCij 
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where |Cy| denotes the area of the cell Cy, |e| the length of the edge e and Wij the average value of the 



unknown on the cell CV 



II 



W(a;, y)dxdy. 



Here, the elementary flux 4>|- across edge e is computed by a Roe type scheme $ e = A+Wjj + A~Wk y i, 
where n = [n x , n y ) is the outward normal to the the edge e, A n = n x A\ + tiyAi and C^; is the 
neighboring cell of Cy sharing the edge e with it. We can rewrite l|5.1|l as 

/- Wij \Ai\Wjj + AjWj+uj - A+Wj-u , l^2|Wi.j + - _ 



Ax 



Ay 



Ax - Ay 

We will further denote Ax — —— and Aw = — r— , the non dimensioned counterpart of the mesh size in 

cAt cAt 

x and y directions. 

5.2 Discretizations of the interface conditions 

In the following we consider the discretization of the interface conditions. In order to do that we will first 
write the semi-discrete system (only in the y direction): 



(5.3) 



where 



pPj + ud x P 3 + V my P 3 + pc 2 d x Uj + pc 2 V py V 3 = fx 
(3U 3 + ud x U, + vD-U 3 + -_d x P, = f 2 



(3V, + ud x Vj 



T> — 



_D py Pj — J3 



c + v^_ c 

— D ,+- 



where we M n = ^ and M t = S are the normal and the tangential Mach number and are the usual 
finite difference operators. When writing the discrete counterpart of the algorithmic there is no need for 
a special treatment in the update step since the boundary conditions do not involve any derivatives. But 
one has to take care of the discretization of the interface conditions in the correction step. For this, we 
proceed as in the proof of proposition [21 but at the semi-discrete level l|5.3f) . First of all, from 1)5. 3[) we 
obtain the discrete counterpart of l|4.13|l 



(5.4) 



d n P = 



:[-Un(P + u T V my ){P) + pc 2 {P + u T D-){U n ) - u n pc 2 V py (U T )} 



and then by replacing it into l|4.11|l we obtain 



(5.5) 



B d is{W) = u n pc 2 [-(/3 + u T D ){U n ) + u n V py {U T ) + u n (u T d y - V my )P] 



which suggests the use of the classical finite difference operator D y for the normal components of the 
velocity and T> py for the tangential one. Nevertheless, the question that remains open is which kind of 
approximation we should use for the pressure term. The most natural one is to cancel the pressure term. 
We will see in the discrete analysis that follows that it may lead to a diverging algorithm. Therefore, we 
will not completely drop this part but consider also the discretization of (u T d y — T> my )P by AyD+ D~ P. 
As we shall see below, it leads to a converging algorithm. We say then that we have a stabilized algorithm. 
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Theoretical discrete convergence results We will proceed to a discrete convergence analysis as in 
DN04I in order to decide which discretization of the boundary conditions is better. We will recall briefly 
the key ingredients of this analysis. We perform a discrete Fourier transform, by looking for the solution 
under the form: 

3 

(5.6) Wij =J2Y1 ««ze^2)M?)AV^y ; (£) 

£ l=l 

where I 2 = — 1. By introducing this expression into the discrete equation we get that for each £, A;(£) 
and Vj(£) have to be the solution of 

M+ Ar + Ay )Vi®=0. 

If we denote by = - — ^ > ~ 1 and by e y (£) = e , we can further solve this systems numerically 

for each wave-number £. By introducing these expressions in the interface conditions, we get the discrete 
convergence rate. 



Figure 1: Convergence rate vs. Fourier number £ for different values of the normal Mach number, no 
stabilization used 



Figure 2: Convergence rate vs. Fourier number £ for different values of the normal Mach number, with 
stabilization term 

In figures [T] and |2| we show for a flow normal to the interface the convergence rate as a function of the 
wavenumber £ along the interface. Two possibilities arc considered: no stabilization and the stabilization 
mentioned above. These figures show the we need of the stabilization term for the high wavenumbers. 

6 Numerical results 

We compare the proposed method with the stabilization term and the classical method defined in [CFS98 
or DLN04 involving interface conditions that are derived naturally from a weak formulation of the 
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M n 


Classical 


New DDM 


Mn 


Classical 


New DDM 


0.001 


67 


18 


0.4 


25 


14 


0.01 


66 


16 


0.5 


20 


14 


0.1 


55 


14 


0.6 


16 


14 


0.2 


41 


16 


0.7 


13 


12 


0.3 


32 


16 


0.8 


15 


14 



Table 1: Subdomain solves counts for different values of M, 



underlying boundary value problem. We present here a set of results of numerical experiments on a 
model problem. We consider a decomposition into different number of subdomains and for a linearization 
around a constant or non-constant flow. The computational domain is given by the rectangle [0 , 4] x [0 , 1] 
with a uniform discretization using 80 x 20 points. The numerical investigation is limited to the solving 
of the linear system resulting from the first implicit time step using a Courant number CFL=100. For 
all tests, the stopping criterion was a reduction of the maximum norm of the error by a factor 10~ 6 . We 
consider first a decomposition into 2 subdomains and a linearization around a constant state and we are 
solving the homogeneous equations satisfied by the error vector at the first time step. In the following, for 
the new algorithm, each iteration counts for 2 as we need to solve twice as much local problems than with 
the classical algorithm. For an easier comparison of the algorithms, the figures shown in the tables are the 
number of subdomains solves. Table ^ summarizes the number of subdomains solves for different values 
of the reference Mach number for the new and the classical algorithm for a normal flow to the interface 
Mt = 0.0. The same results are presented in figure 03 In Table we consider a linearization around a 

Subdomain solves for 2 subdomains, Mt=0.0, eps=1 .Oe-6 

70 



60 
50 
| 40 
30 
20 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Mn 

Figure 3: Subdomain solves counts for different values of M n 

variable state where the tangential velocity is given by Mt{y) = 0.1(1 + cos(?t?/)) and we vary the normal 
Mach number. In figure 01 we linearize the equations around a variable state for a normal flow to the 
interface (M t = 0.0), where the initial normal velocity is gives by M n (y) — 0.5(0.2 + 0.04tanh(y/0.2))). 
The sensitivity to the mesh size is shown in the Table |3| We can see that for the new algorithm the 
growth in the number of iterations is very weak as the mesh is refined, the same property being already 
known for the classical one. The next tests concern a decomposition into 3 subdomains. We first consider 
a linearization around a constant state. Table Ogives the results for different values of the reference Mach 
number for the new and the classical algorithm. The same results are presented in figure [SJ In figure 
|S| we consider for a three subdomains decomposition a linearization around a variable state for a normal 
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14 
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0.3 


20 


14 


0.8 


14 
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Table 2: Subdomain solves counts for different values of M n , M t (y) 



Convergence history: classic vs. new algorithm 




Figure 4: Convergence curves for the classical and the new algorithms 
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59 


16 



Table 3: Subdomain solves counts for different mesh size 



Table 4: Classical vs. the new algorithm 
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Table 5: Subdomain solves counts for different values of M, 
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Iteration number 3 subdomains, Mt^O.O, eps=1 .Oe-6 



Classical Ovr. 
New algorithm 
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Mn 

Figure 5: Convergence curves for the classical and the new algorithm 

flow to the interface. The normal velocity is given by M n (y) — 0.5(0.2 + 0.04tanh(y/0.2))) (the same as 
for the 2 subdomain case). These tests show that the new algorithm is very stable with respect to various 



Convergence history: classic vs. new algorithm 




Figure 6: Convergence curves for the classical and the new algorithm 

parameters such as the mesh size and the Mach number. We see that the convergence in two iterations of 
the continuous algorithm is lost at the discrete level although the subdomain solves are very reasonable. 
Moreover, a stabilization was necessary for the discretization of the interface condition l|4.3l) in order to 
keep the algorithm converging. The optimal discretization of this interface condition is not yet quite well 
known. The comparison with the classical algorithm is favorable for Mach numbers smaller than 0.5 and 
especially very low Mach numbers by a factor of almost 4. 
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7 Conclusion 



In this paper we designed a new domain decomposition for the Euler equations inspired by the idea of 
the Robin- Robin preconditioner |AJN97| applied to the advection-diffusion equation. We used the same 
principle after reducing the system to scalar equations via a Smith factorization. The resulting algorithm 
behaves very well for low Mach numbers, where usually the classical algorithm doesn't give very good 
results. We reduce the number of subdomain solves by almost a factor 4 for linearization around a constant 
and variable state as well. A general theoretical study and more comprehensive numerical tests have to 
be done in order to firmly assess the applicability of the proposed algorithm to large scale computations. 
This work can also be seen as a first step for deriving new domain decomposition methods for the 2D 
or 3D compressible Navier-Stokes equations. Indeed, the derivation of the algorithm which is based on 
the Smith factorization is in fact general and can be applied to arbitrary systems of partial differential 
equations, see |NR| for the Stokes and Oseen equations and [DNR for a general presentation. 
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